Regresión Lineal Simple: Enfoque Machine Learning

Estas notas estan basadas en el libro R for Data Science de Garrett Grolemund y Hadley Wickham

Modelo Lineal: Enfoque de Machine Learning

library(tidyverse)
library(modelr)
library(plotly)
options(na.action = na.warn)

En el contexto del curso se observa el enfoque estadístico para hallar las estimaciones de los parámetros en un modelo lineal. Obtenemos los valores con las siguientes fórmulas:

\(\hat{\beta_0} = \overline{Y} - \hat{\beta_1} \cdot \overline{X}\)

\(\hat{\beta_1} = \frac{\sum{(X_i - \overline{X})(Y_i - \overline{Y})}}{\sum{(X_i - \overline{X})^2}}\)

En este caso vamos a desarrollar otra manera de hallar los valores de estos parámetros, mediante un enfoque de optimización.

datos sim1

El paquete modlr viene con un set de datos de juguete llamado sim1

ggplot(sim1, aes(x, y)) + 
  geom_point() +
  theme_bw()

Se puede ver un patrón fuerte en los datos. Pareciera que el modelo lineal y = a_0 + a_1 * x puede ser apropiado para modelar estos datos

Modelos al azar

Para empezar, generemos aleatoriamente varios modelos lineales para ver qué forma tienen. Para eso, podemos usar geom_abline () que toma una pendiente e intercepto como parámetros.

models <- tibble(
  a1 = runif(250, -20, 40), # Función para tomar valores al azar de una distribución uniforme
  a2 = runif(250, -5, 5) # Función para tomar valores al azar de una distribución uniforme
)

ggplot(sim1, aes(x, y)) + 
  geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +
  geom_point() +
  theme_bw()

A simple vista podemos apreciar que algunos modelos son mejores que otros. Pero necesitamos una forma de cuantificar cuales son los mejores modelos.

Distancias

Una forma de definir mejor es pensar en aquel modelo que minimiza la distancia vertical de la recta (el modelo) con cada punto:

Para eso, elijamos un modelo cualquiera:

\[ y= 7 + 1.5*x\] Y lo representamos gráficamente

dist1 <- sim1 %>% 
  mutate(
    dodge = rep(c(-1, 0, 1) / 20, 10),
    x1 = x + dodge, # para que se vean mejor las distancias, corremos un poquito cada punto sobre el eje x
    pred = 7 + x1 * 1.5
  )

ggplot(dist1, aes(x1, y)) + 
  geom_abline(intercept = 7, slope = 1.5, colour = "grey40") +
  geom_point(colour = "grey40") +
  geom_linerange(aes(ymin = y, ymax = pred), colour = "#3366FF") +
  theme_bw()

La distancia de cada punto a la recta es la diferencia entre lo que predice nuestro modelo y el valor real

Función de Modelo

Para computar la distancia, primero necesitamos una función que permita representar a cualquier modelo lineal simple. Creamos una una función que recibe un vector con los parametros del modelo, y el set de datos, y genera la predicción:

model1 <- function(ordenada, pendiente, data) {
   ordenada + data$x * pendiente
   }
model1(7, 1.5, sim1)

 [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0
[25] 20.5 20.5 20.5 22.0 22.0 22.0

Necesitamos una forma de cuantificar el error de predicción de nuestro modelo en todas las observaciones y resumirlo en única métrica. Esta será la función objetivo que buscaremos minimizar.

Una de las formas de hacerlo es con el error cuadrático medio (ECM)

\[ECM = \frac{\sum_i^n{(\hat{y_i} - y_i)^2}}{n}\]

# Creamos una función que recibe un vector con los parametros del modelo y el set de datos y devuelve el ECM
measure_distance <- function(ordenada, pendiente, data) {
 diff <- data$y - model1(ordenada, pendiente, data)
 mean(diff ^ 2)
   }
measure_distance(7, 1.5, sim1)

[1] 7.103355

Evaluando los modelos aleatorios

Ahora podemos calcular el ECM para todos los modelos del dataframe models. Para eso utilizamos el paquete purrr, para ejecutar varias veces la misma función sobre varios elementos.

MAP

Nosotros tenemos que pasar pasar los valores de a1 y a2 (dos parámetros –> map2), pero como nuestra función toma sólo uno (el vector a), nos armamos una función de ayuda para wrapear a1 y a2

# Calculamos las distancias para los nuevos modelos
models <- models %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, measure_distance, sim1))
models

A continuación, superpongamos los 10 mejores modelos a los datos. Coloreamos los modelos por -dist: esta es una manera fácil de asegurarse de que los mejores modelos (es decir, los que tienen la menor distancia) obtengan los colores más brillantes.

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(data = filter(models, rank(dist) <= 10),
              aes(intercept = a1, slope = a2, colour = -dist)) +
  theme_bw()

NA

También podemos pensar en estos modelos como observaciones y visualizar con un gráfico de dispersión de a1 vsa2, nuevamente coloreado por -dist. Ya no podemos ver directamente cómo se compara el modelo con los datos, pero podemos ver muchos modelos a la vez. Nuevamente, destacamos los 10 mejores modelos, esta vez dibujando círculos rojos debajo de ellos.

ggplot(models, aes(a1, a2)) +
  geom_point(data = filter(models, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist))  +
  theme_bw()

Grid search

En lugar de probar muchos modelos aleatorios, podríamos ser más sistemáticos y generar una cuadrícula de puntos uniformemente espaciada (esto se denomina grid search). Elegimos los parámetros de la grilla aproximadamente mirando dónde estaban los mejores modelos en el gráfico anterior.

# Crear la grilla
grid <- expand.grid(
  a1 = seq(-5, 20, length = 25),
  a2 = seq(1, 3, length = 25)
  ) %>% 
  # Calcular la distancia
  mutate(dist = purrr::map2_dbl(a1, a2, measure_distance, sim1))

grid %>% 
  ggplot(aes(a1, a2)) +
  geom_point(data = filter(grid, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist)) +
  theme_bw()

Cuando superponemos los 10 mejores modelos en los datos originales, todos se ven bastante bien:

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(grid, rank(dist) <= 10)) +
  theme_bw()

Superficie del ECM

Podemos pasar del gráfico de la grilla de puntos a graficar los mismos datos en tres dimensiones. En el plano xy tendremos a ambos parámetros y en el eje z observamos el valor del del error cuadractico medio (ECM).

Notemos que ya no estamos trabajando con la distancia sino que estamos graficando la superficie del ECM como función de ambos parametros.

Por la fórmula del ECM esta superficie es convexa y presenta un mínimo global.

# Matriz para el grafico
rss_matrix <- matrix(models_grid[["dist"]],nrow = length(b1),ncol = length(b1), byrow = TRUE)

# Grafico usando plotly
rss_graph = plot_ly(x=b0, y=b1, z=rss_matrix) %>% add_surface(contours = list(
  z = list(
    show=TRUE,
    usecolormap=TRUE,
    highlightcolor="#ff0000",
    project=list(z=TRUE)
  )
), reversescale=TRUE)  %>%
  layout(
    title = "Superficie del ECM",
    scene = list(
      xaxis = list(title = "a0"),
      yaxis = list(title = "a1"),
      zaxis = list(title = "RSS")
    ))

rss_graph

Óptimo por métodos numéricos

Podríamos ir haciendo la cuadrícula más fina y fina hasta que nos centramos en el mejor modelo. Pero hay mejores formas de abordar ese problema como herramientas de optimización numéricas. Algunas de ellas son:

• Método Nelder-Mead
• Búsqueda Newton-Raphson.
• Gradient Descent

La intuición de este último método es bastante simple: Se elige un punto de partida en la superficie de la función objetivo y se busca la pendiente más inclinada. Luego, desciende por esa pendiente un poco, y se repite una y otra vez, hasta alcanzar la condición de corte.

En R, podemos utilizar la función optim (). Cuenta con los siguientes parámetros:

• par: vector de puntos iniciales. Elegimos el origen por poner cualquier cosa
• fn: función objetivo, y los parámetros que nuestra función necesita (data)

# Adaptamos la función del ECM para usarla con la optimizacion
measure_distance_opt <- function(params, data) {
 diff <- data$y - model1(params[1], params[2], data)
 mean(diff ^ 2)
   }

# Corremos la optimizacion
best <- optim(c(0, 0), measure_distance_opt, data = sim1)

best

$par
[1] 4.222248 2.051204

$value
[1] 4.529154

$counts
function gradient 
      77       NA 

$convergence
[1] 0

$message
NULL

# Obtenemos los mejores parametros
best$par

[1] 4.222248 2.051204

# Graficamos la linea de los mejores parametros
ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(intercept = best$par[1], slope = best$par[2], color="firebrick") +
  theme_bw()

Para llegar a este resultado definimos:

• Tipo de modelo: consideramos que nuestro conjunto de datos (fenómeno) se puede modelar con una regresión lineal simple
• Función de costo: utilizamos el error cuadrático medio
• Método de optimización: la función optim() usa por default el método Nelder-Mead

Función lm()

Al comparar el valor de los parámetros obtenidos mediante el método de optimización con los coeficientes estimados por la función de modelo lineal, observamos que son los mismos.

summary(lm(y~x, data=sim1))

Call:
lm(formula = y ~ x, data = sim1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1469 -1.5197  0.1331  1.4670  4.6516 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.2208     0.8688   4.858 4.09e-05 ***
x             2.0515     0.1400  14.651 1.17e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.203 on 28 degrees of freedom
Multiple R-squared:  0.8846,    Adjusted R-squared:  0.8805 
F-statistic: 214.7 on 1 and 28 DF,  p-value: 1.173e-14

Nota técnica La función lm() llama a la función lm.fit(), la cual utiliza factorización o descomposición QR para resolver las ecuaciones del modelo lineal.